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SINGULAR QUADRATURE RULES AND FAST CONVOLUTIONS 
FOR FOURIER SPECTRAL METHODS 

JAE-SEOK HUH AND GEORGE FANN 

Abstract. We present a generic scheme to construct corrected trapezoidal 
rules with spectral accuracy for integral operators with weakly singular kernels 
in arbitrary dimensions. We assume that the kernel factorization of the form, 
K = aij> + K with smooth a and K, is available so that the operations on the 
smooth factors can be performed accurately on the basis of standard Fourier 
spectral methods. To achieve high precision results, our approach utilizes the 
exact evaluation of the Fourier coefficients of the radial singularity <j>, which 
can be obtained in arbitrary dimensions by the singularity isolation/truncation 
described in this article. We provide a complete set of formulae for singularities 
of the type: log(r) and r~ Convergence analysis shows that the constructed 
quadrature rules exhibit almost identical rate of convergence to the trape- 
zoidal rule applied for non-singular integrands. Especially, for smooth data, 
the corrected trapezoidal rules converge supcr-algcbraically. 



1. Introduction 
We consider the problem of evaluating the integral, 

(1-1) ('/)(*)= / K(x,y)f(y)dy, 

J D 

where the kernel K may have a point singularity at x = y. We assume that / 
is compactly supported on D and can be represented satisfactorily by a Fourier 
series. When the kernel is also smooth (or sufficiently regular for the purpose), 
the usual trapezoidal rule with the integrand (K ■ /) evaluated on a uniform grid 
is a classical rule of thumb for the construction of a spectral scheme involving 
the integral operator (1.1). Additionally, many kernels of interest are given by 
K(x,y) = K(x — y). Then, the Fourier representation of the kernel and the data 
enables the fast evaluation of the multiplier operator I via the FFT. However, the 
singularities in most kernels of interest render the trapezoidal rule unapplicable, 
and we need a corrected one. 

We assume that the kernel can be factored to the form, 

(1.2) K(x,y) = a{x,y)4>{r) + K{x,y) 

where a and K are smooth functions and r = \\y — x\\ is the Euclidean distance in 
W 11 . We assume that the singularity is known and carried entirely by the radial 
function <j>. For the derivation of the quadrature weights, a and K may depend 
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explicitly on the target point x, in which case the weights should be constructed for 
each x. The main subject of this paper is the construction of corrected trapezoidal 
rules which exhibit the same spectral accuracy as the usual trapezoidal rule applied 
for smooth kernels. 

The constructed quadrature rules can be written as 

(1.3) (If)(x) ~ h J2 + h E K(x,y e )f(y e ), 

yeeB R yi=£x 

where the contribution of the singularity is reproduced by the correction weights 
W£ in a neighborhood Br of x. The presented scheme can also be viewed as a 
regularization of the kernel such that the replacement of 4> with the regularized 4> 
does not impair the accuracy of the regular trapezoidal rule with a smooth kernel 
at the given sampling frequency. In what follows are a few comments regarding 
the key features of our method and the relationship to other quadrature rules and 
applications available in the literature. 

(a) One of the main advantages of our approach is the high precision of the result- 
ing quadrature rules. The construction does not involve any other non-trivial 
quadrature rule or a solve of a linear system both of which typically limit the 
achievable order of accuracy. In our approach, the FFT is the only numeri- 
cal algorithm used and the accumulation of the rounding error of the FFT is 
only 0(e mac h log N). To achieve the high precision, the Fourier coefficients of 
the truncated singularity (j> are evaluated up to the machine precision. One of 
the major contributions of this paper is the complete recipe we present for the 
logarithmic and the power-law singularities. 

(b) The second advantage is its applicability for any dimension. The key functions 
required for the construction can obtained by recurrence relations presented in 
main body of the paper. Except for the r~ v singularities with a non-integer v 
(which is less significant in applications than the integer cases, we think), the 
construction does not require an implementation of a new special function. 

(c) We present the factored forms of the Hclmholtz kernels in arbitrary dimensions 
so that the corrected quadrature rules can be readily obtained from them. The 
Hclmholtz kernels are of great importance in applications. We present numer- 
ical examples which demonstrate the advantage of the presented quadrature 
rules in the construction of high order integral operators especially for oscillat- 
ing kernels with high wavenumbers. 

(d) We do not pursue in this study the end-point corrected trapezoidal rules for 
non-periodic data, which can be found in [2, 14, 19]. In principle, an end-point 
correction is equivalent to an accurate estimation of derivatives at and near the 
end-points, which can be done by an over-sampling or by a certain type of data 
extension. For the latter, a simple extrapolation results in a terrible oscillation, 
which limits the achievable order of accuracy. Hence, most successful high-order 
methods involves a solve of a least-squares problem. An accurate and stable 
extension is possible also for the trigonometric basis functions (cf. [13]). We 
can represent a smooth non-periodic function by a smooth compactly supported 
interior and a near boundary part. In this paper, we focus on the former case 
and the correction for the near boundary part will be considered in a separate 
study. The quadrature rules in [2, 14, 19] are for one-dimensional singular 
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integrals and the multidimensional extension is not so straightforward since 
a weak singularity in a higher dimension can not always be represented by 
a product of one-dimensional weak singularities. Multi-dimensional singular 
quadrature rules of relatively low orders can be found in [21], which are applied 
for non-oscillatory kernels. 

(e) When / is periodic and the domain of the convolution (1.1) is unbounded, the 
kernel can be factored as K = <p K + (1 — <p)K by a smooth cut-off function (p. 
The convolution can be recast to 

(1-4) // = J K ■ (p/) +/((!- f)K) ■ /, 

where the first integral involves the singular kernel and the compactificd ((pf), 
and the second integral the convolution of the periodic / with the regularized 
kernel (1 — <p)K where the trapezoidal rule can be applied without any cor- 
rection. The Ewald summation (cf. [10]) focuses on the design of the smooth 
kernel so that the Fourier coefficients of the regularized kernel decays rapidly. 
Classical applications of the Ewald summation are for point sources, hence, 
the first integral is merely a summation. Our focus is more on the accurate 
evaluation of the first integral for continuous /. 

(f) The data / or their Fourier coefficients do not necessarily have to be given on a 
uniform grid as long as equivalcntly fast and accurate algorithms are provided. 
For the extension of the quadrature rules to non-uniformly sampled data, one 
can refer to [11, 12]. For sparse data as well as rapidly decreasing kernels, the 
use of the partial FFT (cf. [3]) can also be utilized. 

(g) The kernel factorization of the form (1.2) can be found in numerous articles. 
Probably, the most well-known examples are the Nystrom methods for bound- 
ary integral equations such as [9, §3.5]. In [7, 22], combined with partition-of- 
unity boundary decomposition, Nystrom-type schemes are presented for bound- 
ary integral equations in M 3 . Spectral schemes in [7, 22] commonly remove the 
1/r-singularity by a polar change of variables, and evaluate each of the resulting 
one-dimensional integrals by the trapezoidal rule using resampled data without 
building quadrature weights explicitly. 

(h) In order to achieve a proper spectral convergence, the kernel should be cor- 
rectly factored so that a and K arc both smooth. Regarding this issue, we 
should keep the fact in mind that r 2 is smooth but r is not even in C 1 . The 
following factorization might look attractive since the singularity is isolated 
independently of the wavenumber. 

,„ e tkr 1 f cos(kr) - 1 sin(/cr)\ 

(1.5) = h — h i — — - . 

4irr Airr \ Airr Airr J 

However, the above form should be avoided since (cos(Ar) — l)/r is not 
smooth. A correct form of factorization is 

e lkr cos(fcr) .sin(fcr) 
Airr Attt 4-rrr 

where a(r) — cos(fcr)/(47r), 4>{r) = 1/r, and K{r) = i sm(kr) / (4irr) . Readers 
will soon notice that all the well-behaved functions appearing in this paper have 
power series expansions containing only even powers of r. 
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This paper is organized as follows: in §2, we begin with the detailed derivation 
of the quadrature rules, which will be followed by a complete set of formulae for 
the evaluation of the Fourier coefficients of the radially-truncated singularities of 
the type: log(r) and r~ v . Then, we present the factored forms for the Hclmholtz 
kernels in arbitrary dimensions. In §3, we discuss the order of accuracy of the quad- 
rature rules related to the regularity of the data. Finally, the results of numerical 
experiments will be presented in §4. 

2. Construction of Quadrature Weights 

2.1. Notations, grids, and DFTs. Let U = [— 1, l] m be the computational do- 
main in W 1 . We consider a uniform grid on M m with the spacing 1/Nj in the jth 
dimension. Thus, for any multi-index I = (£\, . . . , £ m ) € Z m , the corresponding 
grid point yi = (yg 1 , . . . ,ye m ) is given by ye- — tj/Nj. We assume that the physical 
domain is given as an affine image of U (without translation for the sake of sim- 
plicity), that is, a parallelotope in M. m . We denote by \ this non-degenerate linear 
mapping U i— > then the related metric tensor is constant and is given by x T X- 
We use the symbol r for the Euclidean distance in the physical domain. That is, 

(2-1) r(x,y) = \\ X (y-x)\\. 

With the target x £ U fixed and r viewed as a function of y e U, we view (1.1) 
as an integral on U with the Jacobian determinant \x\ = \J dct(x T x) multiplied 
afterward. The function r is not smooth at the point of singularity due to the square 
root but r 2 is a smooth quadratic polynomial. Note also that the singularity <p is 
not radial in U (but is radial in xU)- 

Let N ee rjJLi Nj and define l N C Z m by 

(2.2) l N ee {(4, ...,l m ) \ -Nj< tj < Nj - 1}. 

Let F be a periodic function on U . Application of the trapezoidal rule to the 
(periodic) Fourier transform 

(2-3) %) = i / F(y)e-™ k -ydy 

results in the formula for the DFT 

(2-4) ^ = ^E J? («) e_i "*" W - 

The truncation applied to the inverse Fourier transform (of delta functions on Z m ) 

(2.5) F(y)= J2 F(k)e M -y 

feez™ 

defines the inverse DFT by 

(2.6) Fi — Fk e i7rk ' ye . 

fceijv 

Note that the Fourier coefficients F k approximated by the DFT contains error 
due to the aliasing. Suppose we are given exact samples, i.e. Fg — F(ye). The 
interpolation F of the samples defined by 

(2.7) F(y) = Fk e™ k ' y where F k is the DFT of F t 

fceijv 



SINGULAR QUADRATURE RULES AND FAST CONVOLUTIONS 



5 



involves two sources of error: (1) the aliasing contained in Fk and (2) the truncation 
error of the inverse DFT. The decay characteristics of the exact Fourier coefficients 
F(k) takes a central role in the estimation of the error. 

2.2. The support of the data. We assume that the data (or source) / in (1.1) 
viewed as a function on the computational domain is supported on [0, l] m . Then, 
/ can be extended on U by padding zeroes. The rationale behind this is quite 
obvious; since we take x + U as the support of the kernel K(x, •) (or equivalcntly, 
the domain of integral) for the target x, if / is not extended, the convolution will 
include the effect of the fictitious portion of the periodized source. Similarly, if x 
is not in [0, l] m , the result (If)(x) will include the effect of the fictitious source, 
hence, the portion on U \ [0, l] m of (If )(x) obtained on U should be discarded. 

For x ^ [0, l] m , the integrand becomes smooth, hence, the usual trapezoidal rule 
serves our purpose well without any correction. The resulting summation can be 
accelerated by any fast multipole (or an equivalent) method designed for the kernel. 
In this paper, we focus ourselves on the near-field solution for x G [0, l] m . 

On this setting, the translation f(y — x) by x £ [0, l] m results in the translated 
source vanishing on the boundary of U and satisfying all the assumptions. There- 
fore, in this section, without loss of generality, we assume that the target point x 
is at the origin. Every function is viewed as a function of the source point y only 
and the symbol x will be omitted. Thus, we recast (1.1) to the integral, 

(2-8) 1(f) = \X\ I K{y)f{y)dy. 

Ju 

2.3. Localization of the singularity. Let r(y) = \\xu\\- With a slight abuse of 
notation, we denote by Br an ellipsoid in U given by Br = { r(y) < R} for some 
R > 0, where R is chosen such that the image X-Br (a ball of radius R in the 
physical domain) is contained in \ U. For the accuracy, the best choice of R is 
m&x y€ gu r(y) so that Br contains as many grid points as possible. 

Let if be any smooth even function on R such that (1) <p(0) = 1, (2) ip = on K\ 
(-R, R), and (3) ip c = (1 — ip) vanishes smoothly at the origin. Then, (<j>(r) (p c (r)) is 
also a smooth function and vanishes at the origin. For all the numerical experiments 
in this paper, we utilized <p(r) — p\(r/R) where <pi is the sigmoidal function, 



(2.9) 




- 2/|t| /(i-|t|) 2 |f| < i 
\t\ > 1 



Then, we isolate the singularity by using the identity 
(2.10) 0(r) = ^(r)p(r) + q)(r)p c (r) 

and rewrite (2.8) as 

(2.H) I(.f) = \x\[ Hr)F(y)dy+\ X \ [ G(y)dy 

Jb r Ju 

where F and G are regular periodic functions given by 

(2.12) F(y) = <p(r)a(y)f(y) 

(2.13) G(y) = V c (r) 0(r) a(y) f(y) + K(y) f(y). 
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The second integral without the singularity can be treated well by the usual 
trapezoidal rule; that is, 

(2.14) f G{y) dy~^Yl a ^ + /(*«)■ 

For the first integral with the singularity, we utilize the interpolation (2.7) to obtain 



(2.15) \ X \f <f>(r)F(y)dy~\ X \Y, ^ / 4>(r) e™^ dy 

(2-i6) = = E nvt) E / *m e<OTfe ^ rf y 

cel.. ten., \ •'Br 



(2-17) ^ W E E 



IXI 



(2-18) =^^F{y e )4, t , 
where we defined 

(2.19) 0(fc) = _L/" ttfe-^-vdy and & = £ 0(fe). 

We can change the sign in the exponential function arbitrarily since r(—y) = r(y). 
Notice that (f)(k) are the exact Fourier coefficients of 4>(r) which is truncated to 
zero on the exterior of Br and viewed as a periodic function on U. Or equivalcntly, 
we can consider cfr as the (non-periodic) Fourier transform (scaled by 2~ m ) of <j> 
truncated on M. m \ Br. By definition, (j>i are the inverse DFT of the finite samples 
4>(k) in the frequency domain. Consider the interpolation 

(2.20) 4>{y)= ]T me^-y. 

Then, 4>i = 4>{yt) and (2.18) is simply the trapezoidal rule applied to the product of 
the two functions, <f> and F. Thus, the procedure is equivalent to the regularization 
of the singular <f> such that the regularized kernel <j> results in the exact integral by 
the trapezoidal rule if F can be exactly represented on the given grid. 

Care should be taken not to confuse <j> with the Fourier transform of the original 
non-truncated (p. To be more precise, we should use a notation like <pn, but we 
choose the notational simplicity. Typically, is a slowly decaying function with the 
point singularity at the origin, hence, without the truncation, its Fourier transform 
possesses the same nature in the frequency domain. By truncating in the space, 
we regularize the Fourier transform to a smooth (but still slowly decaying) <j>R. 
By truncating in the frequency domain, that is, by sampling only up to the given 
sampling frequency, we obtain (f> regularized in the space. 

As one can notice from the above derivation, our construction requires the exact 
evaluation of <f>. In §2.5, we present a detailed discussion on the nature of (j) as well 
as the explicit formulae for the logarithmic and the power-law singularities. 
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(2.25) w e 



2.4. The corrected trapezoidal rule. Merging (2.14) and (2.18), the quadrature 
rule can be written as 

(2-21) /(/) ~ M £ (a(y e ) (<p(r t ) fa + ^( n ) fan)) + K(y e )) f(y e ). 

e&N 

For £ = 0, since ip(0) = 1 and </5 c (0) = 0, the weight becomes 

(2.22) a(0) fa + K(0). 
For £ 7^ 0, we can rewrite 

(2.23) <p(r t ) fa + <y5 c (r £ ) fan) = (fa - fan)) p(ry) + 0(r/), 

Thus, we present the quadrature rule by the sum of the trapezoidal rule and the 
correction rule; that is, 

(2-24) /(/) ~ M £ A- (|tt) + M £ W£ 

where the correction weights W£ are given by 

r a(0) fa + K(0) £ = 

a{yi) (4>i ~ <t>(rtf) <p(rt) £^®' 

The separated representation (2.24,2.25) is the final form of the quadrature rule 
we present in this paper. The quadrature weights (2.25) is quite self-explanatory; 
at the point of singularity, <j> has been represented by the equivalent finite weight 4> 
which accompanies the balancing neighbors (fa — <j>(ri))ip(re) to achieve the desired 
spectral accuracy. The advantage of using the separated form (2.24,2.25) over the 
primitive one (2.21) is obvious; oftentimes, the smooth remainder term K has a 
quite complicated form. Since the separated form requires only the limiting value 
of K at the origin, the implementation can be simpler and more readable. 

There must be a few careful readers concerning the numerical soundness of the 
expression (fa — fare)). With finite precision arithmetics, it sometimes comes to a 
catastrophic end to take the difference of two potentially large values; the cancel- 
lation error can corrupt the result. However, on uniform grids where grid points 
are not clustered, the issue is not so significant. Moreover, for spectral methods 
like ours, when applied to quite smooth data, the solution often converges to the 
desired precision before the grid spacing becomes small enough to bring such an 
issue to the surface. In all the numerical experiments we conducted including all 
the examples in this paper, we utilized the separated representation but have never 
experience a trouble. Actually, the classical Nystrom methods frequently use the 
expression K = K — a (f> to avoid the explicit evaluation of a complicated K (see, 
for example, [9, §3.5]). For rare cases when the issue becomes significant, one can 
utilize the original form (2.21) for grid points close to the singularity. 

We can also bring the convolution (2.24) to the frequency domain. In this sense, 
our construction is equivalent to the representation of the Fourier transform of K 
truncated onR m \[/ by 

(2.26) K~a*fa*if + [a*(j>*(f c + K 
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where the terms in the parentheses decays rapidly. The key idea is that the appli- 
cation of the radial cut-off function ip enables us to truncate <p outside of Br and 
we can evaluate the Fourier transform (which is now the one-dimensional Hankcl 
transform) of the truncated <f> exactly. Numerically, we evaluate K by the DFT of 
{wo, wg + K(ye)}, where we are obtained by the inverse DFT of 4>{k). 

Now, the true meaning of the grid we used so far has become clearer. It is not 
necessary for the data / to be given on that (or on any) grid. Regarding /, the 
only information we need is the required minimum sampling frequency. The grid we 
have used so far is in principle to perform the DFT of K and the convolutions (by 
multiplications in the spatial domain) in (2.26), hence, we call it more specifically 
the construction grid. 

(A) A recommended sampling frequency for the construction (which is greater or 
equal to the given sampling frequency of the data) is such that the interpolation 
errors of the smooth functions a, ip, and K are comparable to that of the data. 
Otherwise, the quadrature weights will still work but the error in the weights 
will be the dominating factor (cf. §4.2 and §4.2). 

(B) For the inverse DFT to obtain <f>, we do not need to use the entire construction 
grid. Only a subset of the grid containing Br is sufficient. The use of the 
subset can be very useful if the kernel exhibits exponential decay like the 
Helmholtz kernel with a complex wavenumber k with Im(fc) > 1. Such an 
exponentially decaying kernel requires a very high sampling frequency for the 
construction (cf. §4.2). However, due to the rapid decay, the kernel outside of 
a certain small Br is practically zero. Hence, we can compute (we + K(ye)) 
only within the small subset of the conceptually huge construction grid. The 
DFT on the entire grid can be performed by padding zeros and the unnecessary 
high frequency terms can be discarded, both of which can be done in a single 
efficient procedure without using huge temporary memory by (so called) the 
partial FFT. 

2.5. The Fourier transform 4>. The construction of the correction weights re- 
quires the exact evaluation of the Fourier transform 4> of the truncated singularity. 
We begin with the definition of <fi, 

(2-27) m = ^[ 4>(r)e-™ k -ydy. 

1 J Br 

Note that </> is not a radial function in the computational domain U and Br is not 
a ball in U. First, we perform the change of variables rj = R~ 1 xy back to a scaled 
physical domain. Utilizing k ■ y — k' ■ rj with k' — Rx~ T k, we obtain 

Dm r 

(2.28) ^ = / ^(RHDe'^'^drj. 

2 Ixl y||„||<i 

The function 0(i?||-||) is radial and supported on the unit ball. Hence, the integral 
on the right side is a radial function of k! and is given by the following Hankel 
transform on one-dimension. 

,2 - 29 » m= ^\HrL *Uv '-"ww* 
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where J M is the Bessel function of the first kind and 

(2.30) Pk =irR\\ X - T k\\. 
Then, we can rewrite (2.29) as 

(2.31) 0(p) = - , ,V — tt— / 01 I mr-^fM, 
V ^ y " 2™|x| T(m/2 + l) 7 v \p J mW 

where the Bessel function of degree (to — 2)/2 has been replaced with a better 
behaved function A m defined by 

(2-32) 

Then, A m (0) — 1 and \A m \ < 1. The series representation of A m can be obtained 
from that of J( m _2)/2- 

,™ „ ^ ^ (-i) £ r(TO/2) /r 21 

(2 ' 33) ^)-E , !r (, + TO/ 2) (2- 

Thus, ^4 m is an entire function on C and is even on R. The asymptotic behavior of 
A m , which can be obtained from that of J( m -2)/2, governs the decay characteristics 
of <f>. For t > to 2 , 

A (A r(m/2) cos(t - (to - 1)tt/4) 

Note that A m includes a Bessel function with an integer index for an even to 
and with a half-integer index for an odd to, the latter of which can be expressed by 
a finite series of trigonometric functions. Although the series representation may 
look complicated, A m actually consists of familiar functions. For the first 4 indices, 

(2.35) A 1 (t)=cos(t) 

(2.36) A 2 (t) = J (t) 

(2.37) A 3 (t) = sm(t)/t = sinc(t) 

(2.38) A 4 (t) = 2J 1 (t)/t. 

For indices greater than 4, closed-forms can be obtained by utilizing the three-term 
recurrence relation, 

(2.39) A m+4 (t) = m{m t2 + 2) (A m+2 (t) - A m (t)) , 

which can be easily verified from the series form. Thus, A m involves only J and 
Ji for an even to, and only cosine and sine for an odd to, which does not introduce 
any implementation issue. 
Another useful identity is 

(2.40) it m A m+2 (t))' = m t m - x A m (t), 

which results in another expression for (2.31); suppose <j> is C 1 on (0, R] and 
(4>(t) t m ~ e ) — > as t — > for some e > 0. Applying the integration-by-part, (2.31) 
becomes 

rz m Dm f pi -< 

(2 ' 41) ^ {P) = 2"| X | r(m/2 + 1) \<t>{R) A m + 2{p)- j Q <t>i{Rt)t m - l A m+2 {pt)dty 
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where (f>i(t) — t(j)'(t) and the change of variables t/p ^ t is applied for the second 
integral. Note that ^(t) = 1 for <j)(t) = log(i) and ^(t) = -v<f>{t) for <j>{t) = t~ v , 
which provides us with useful recurrence relations for the integral term. The formula 

(2.41) also reveals the fact that the asymptotic behavior of <p is determined by that 
of A m+2 (not A m ), which is summarized in the following lemma. 

Lemma 2.1. Let the grid parameters, R and \, be fixed. Suppose <f> is C 1 on (0,R] 
and (<f>(t) t m - £ ) -> as t -> for some e > 0. Then, for k G Z m , 

(2.42) |0(fc)|=o(||fc||-( m+1 )/ 2 ). 
Proof. First, consider the integral in (2.41). 



(2.43) 



l 

m— 1 



MRt)t m -'A m+2 (pt)dt 







< \\Mt)t m \\oo f 1 \Am+2{ P t)\ ^ 



R m In t 



where we denote by || - ||oo the £°°-norm on [0, R}. From (2.34), there is a constant 
C > such that \A m+1 (pt)\ < C(pt)-^ m+1 ^ 2 . Define a by C = (pa)( ro+1 )/ 2 . 
Then, a < 1 if p > C 2 ^ m+1 \ Since \A m+2 \ < 1, 

1 1^+2(^)1 df < ^ a ^_i dt + c p - (m+ i)/2 £ te - (ro+ i)/2-i ^ 

„e „— (m+l)/2 . . 

(2.44) <8- + — f 1 - a-^+D/ 2 ) 

v ; ~ e e- (m + l)/2 V J 

Cp -( m +l)/2 , S 1 



e-(m+l)/2 V e e-(m + l)/2 / 

where <5 is an arbitrary constant > 1. When e > (m + l)/2, a £ = O (p _ ( m+1 )/ 2 ) 
since a = O (p -1 )- When e < (m + l)/2, we can choose S — e/(e — (m + l)/2) > 1 
so that the second term vanishes. Thus, J* t e - 1 \A m+2 (pt)\ dt = O ( p -(™+i)/2) . 
Combined (2.34) for the first term in (2.41), the above result shows that \<j>(p)\ = 
O (p~( m+1 ^ 2 ). Since pk < 7ri?||x _T ||||fc|| where ||x~ T | is the operator norm of 
X~ T , the same relation holds for ||fc|| also. □ 

In what follows, we present the formulae of <j> for <j>(r) = log(r) and for <p(r) = r~ v 
with v < m, which enables us to evaluate them up to the machine precision. We 
begin with the logarithmic singularity. 

2.5.1. 4>{r) — log (r). From the integration-by-parts formula (2.41), we obtain 

(2-45) &(p) = 2ro J r(m/2 + 1) {lo g (i?)^ + 2(p) - L m (p)) 

where 

(2.46) L m (p)= f t™- 1 A m+2 {pt)dt. 

Jo 

For m = 1 and 2, 



(2.47) L 1 (p)= [ 1 MPll dt= m 

Jo Pt P 

(2.48) L 2 (p) = / 

Jo 



1 2 Ji(pt) ^ = 2(1 - J (p)) 
p P 2 
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For m > 3, utilize the recurrence relations (2.39) to obtain 

(2.49) L m+2 (p) = m(m 2 +2) ( t m ~ l (A m+2 (pt) - A m (pt))dt 

P Jo 

where the second term can be integrated by using (2.40), which results in the 
recurrence relation for L m , 

(2.50) L m+2 (p) = i^-t^l | m i ro (p) - A m+2 (p)|. 

With the formulae for L\ and L 2 , the above recurrence relation enables us to con- 
struct an explicit formula of </>^ for any to. In order to avoid possible cancellation 
error of the recurrence formula for small p, we need to evaluate the series represen- 
tation, 

(9V\ T r.^-V (-l)T(m/2 + l) fp\ 2£ 

which is also more efficient for small p. 

2.5.2. 4>(r) — r~ v (v < m). First, notice that r~ v with any v < to can be 
factored as r~ v = r 2n ■ r~ e with a non-negative integer n and e € [to — 2, to). Since 
r 2 is smooth, the smooth factor r 2n can be included in a. Interestingly, due to the 
factor r 2n in a, the correction weight w at origin vanishes for r~ v with v < m — 2. 
For such weaker singularities, the correction is made by the correction weights (with 
relatively smaller magnitude) near the location of the singularity. 

Therefore, from now on, we presume v £ [to — 2, to). It is convenient for the 
presentation to use a new notation, p = m — v, where p e (0,2]. Then, (2.31) is 
written in terms of p as, 

where 

(2.53) M<£\p)= ( mt^ 1 A m {pt)dt. 

Jo 

The series representation is given by 

which can be evaluated efficiently for small p avoiding the cancellation error of the 
recurrence relation. 

We can derive a recurrence relation for M m like the one for L m in the previous 
section. By integrating by parts and utilizing tA' m+2 (t) — m{A m {t) — A m+2 (t)), 

pM^l 2 (p) = (m + 2)A m+2 (p)-(m + 2) pi? A' m+2 (pt) dt 

(2.55) Jo 

= (to + 2) A m+2 (p) - (to + 2) MM(p) + mM^ +2 (p). 

Hence, 

(2.56) (m-p)M^Up) = (m + 2){A m+2 (p)- M#>(p)). 



21 



12 



JAE-SEOK HUH AND GEORGE FANN 



And, the above equation implies that Mm (p) = A m+ 2(p). Since /i e (0,2], there 
are only two cases with \i = m: it = m = 1 and p, = m — 2. Except those two 
cases, we can utilize the recurrence relation, 

(2-57) M<# 2 (p) = ^±l^ m+2 (p) - mM(p)). 

Among /x € (0, 2], the two integer cases are of prime interest; (1) r 2 ~ m (p = 2) 
is the principal singularity of the Helmholtz kernel in M. m (in even dimensions with 
m > 4, the Helmholtz kernel contains an additional logarithmic singularity, which 
can be treated by <f) m in the previous section). Interestingly, this (probably) most 
important class of singularities has the simplest description; M m ' can be written 
explicitly without using the recurrence relation. (2) Singularities with p = 1 arise 
when the Helmholtz kernel in M m+1 is acting on m-dimensional flat boundary. For 
m = 1, the domain of integral need not be a flat manifold (see an example in §4.4 
for the application of the quadrature rule on curves) . We are studying the extension 
of our method for higher dimensional general (non-flat) manifolds. 

(I) ii = 2 {v = m- 2). 

(2.58) MW(p) = [ mtA m (pt)dt 

Jo 

(2) (2) 

The evaluation of M{ ' is straightforward. The formula for M 2 ' is the result 
of (2.56) with fi = m. For m > 2, it is not difficult to obtain the formula 
from the series representation of A m . 



(2.59) 




cos(p) — 1 


— sinc(p) 


P 1 


(2.60) 


Mf 


= Ai(p) 




(2.61) 


M$ 


m(m — 2) 
= P 2 


(l-A m _ 2 (p)j form>2 


(II) /x = 1 


(y — m 


-1). 




(2.62) 




M«(p) = 


= / mA m (pt)dt 
Jo 



Formulae for and requires elementary calculus only. To the best 

of our knowledge, there is not a simple representation of by well-known 
functions. Hence, we treat the following integral form of as the defini- 
tion of a special function (see [1, p. 480] for the properties of the integral). 

(2.63) M[ 1 \p) = A 3 (p) 

(2.64) M 2 (1) (p) = 2/ J Q (pt)dt=- J (t)dt 

Jo P Jo 

(2.65) M^pH 38 ^ 

One may implement his/her own version of from the series representa- 
tion J and the asymptotic expansion of the integral, or can simply use an 
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implementation of f£ J (t) dt in Algorithm 757 (MISCFUN) of ACM Trans- 
actions on Mathematical Software (TOMS) (cf. [17]). For m > 1, Af£} 2 can 



be obtained from the recurrence relation, 
(2-66) Mj» 2 {p) = (V«(p) - A m+2 {p)) ■ 

(III) \x G (0, 1) U (1, 2). 



The recurrence relation (2.57) can be applied to obtain f° r &n Y m> 1. 

Hence, we only need to consider two initial cases m = 1,2. 

(a) m = 1. 



(2.67) 



K (P) = / t^- l cos{pt)dt 



Ci(p,p) 



o 

2_ (-l)V* 



£ (2£+l)\(2£ + p)- 



The function Ci(/x, p) = / Q p t^ 1 cos(t) dt is known as the generalized 
cosine integral, which is related to the lower incomplete gamma function 
with pure imaginary argument, 



(2.68) Ci(p,p)=Re((-i)"7(/MP))- 
(b) m = 2. 

(2.69) M 2 W (p) = f 1 ^- 1 J (pt)dt = ^ { ~ 1)e{p/2) 

Jo 



00 c iW„/o"\2^ 
£(fl)2(2* + / x)' 



To the best our knowledge, there is no available/reliable implementation of 
either M 1 (Ai) or (or any special function which can be used to compute 

them) in the public domain. We can follow the standard implementation 
procedure of special functions - the partial sums of the above power series' 
for small p and the asymptotic expansions for large p. For , the error of 

the 18-term partial sum is less than 1CP 16 on < p < 2ir. For 20 terms 

are enough for the error less than 10~ 16 on < p < 7.01558666981561875 
(the upper limit is the second positive zero of J\). For p > 14it for 
and p > 44.7593189976528217 (the 14th positive zero of Ji), asymptotic 
expansions described in Appendix A and B produce 16-digit accurate results. 
For p in the intermediate range, we divide the domain as U^ =1 [a„, a„ +1 ] 
where a„ = 27m for and a n is the (2n)th positive zero of J\. When 

p G [a n , a n+ i] for some n — 1, ... ,6, 



(2.70) MW(p) = MW(a n ) + l /V 

P J a n 



9 ^-iJcos(t)di m=l 
J (t)dt m = 2 



where Mm\a n ) can be precomputed and reused, and the integrals on [a n , p] 
are evaluated at each time by a quadrature rule. Since the integrands are 
very smooth and oscillate less than one cycle in the interval, any high or- 
der numerical quadrature (such as Clenshaw-Curtis) with a small number of 
samples can compute the result up to the machine precision. 



14 



JAE-SEOK HUH AND GEORGE FANN 



Thus, we have presented our scheme for the construction of corrected quadrature 
weights and the required formulae for the Fourier transforms of logarithmic and 
power-law singularities. We conclude this section by presenting the factored form of 
the Helmholtz kernel in arbitrary dimension, for which the results we have developed 
so far turn out to be well-suited. 

2.6. Helmholtz kernels. Denote by K k (r) the Helmholtz kernel in K™ with com- 
plex wavenumber k such that Im(k) > 0. The domain of the integral (i.e. the 
domain of the convolution) is not necessarily n-dimensional. If the convolution is 
performed on R m , n needs only to satisfy n < m + 1. Let z = kr and v = (n — 2)/2, 
then the Helmholtz kernels are given by 

* V „«>,., 1 1 k \"~ 2 < Y -M , , M*) 



where is the Hankel function of the first kind. More familiar forms in the first 
three dimensions are 

(2-72) JffW = ^, if 2 fc (r) = ^(fcr), K*{r) = ^. 

In the limiting case with k = 0, K® is simply a constant multiple of the logarithmic 
or the power-law singularity; 

(2.73) # 2 °(r) = -^log(r), 

A»=%^ for„>3. 

Recall that J 1/ (z)/(z/2) v = A n (z)/T(n/2) with A n defined and extensively used 
in previous sections. Like A n , the imaginary part of K k is analytic, and has the 
limiting value 

n-2 

ki 



(2.74) limlm(i^(r)) - ... 

Thus, the singularity is carried entirely by the real part. In principle, for any n, 
Kc(K^) contains the same type of singularity as K® given above. However, for even 
n > 3, the kernel contains an additional logarithmic singularity. 

2.6.1. Odd n. Since v = (n — 2)/2 is an half-integer, we can utilize the identity, 
Y v (z) = (-1) M J-u(z), to obtain 



(2 ' 75) {z/2Y ( ^ (z/2)-" (z/2)^ 

(2 76) - -f-nrn/21 { lY'* M-njkr) _1_ 

(ZJb) 1 L) [kj r(2-n/2)r»" 

Therefore, the factored form is given by 

(2-77) K k n {r) = °0+K k n {r) 



2 ' 
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with smooth functions, 

(_l)\n/2] 

(2.78) a k (r) = a k (0)A i _ n (kr) where a k (0) = [ \ 

4r(2 — n/2) ^tt 

(2.79) K k (r)=K k (0)A n (kr) where tf*(0) = -^j^ \^=J . 
It is convenient to use the following recurrence relations. 

(2.80) a k n+2 (Q)= n ^a k n {Q) and K k +2 (0) = ^- K k (0). 

Recall that, in the separated form of the quadrature rule (2.24), we do not need to 
evaluate K k explicitly; only the limiting value at the origin is required. However, 
the values of at grid points are still needed. The following table shows explicit 
formulae and values for the first few dimensions. 



n 


0(r) 


«£(r) 


a k (0) 


^(0) 


1 


r 


— sinc(fcr)/2 


-1/2 


i/(2fc) 


3 


r- 1 


cos(fcr)/(47r) 


1/(4tt) 


zfc/(47r) 


5 


r" 3 


( cos(fcr) + kr sin(/cr)) / (87r 2 ) 


1/(8tt 2 ) 


ifc 3 /(247r 2 ) 


7 


r- 5 ((3- 


(kr) 2 ) cos(fcr) + 'Skr sin(fcr)) /(16tt 3 ) 


3/(16^ 3 ) 


ifc 5 /(2407r 3 ) 



2.6.2. Even n. The Bessel function of the second kind with an integer index v = 
(n — 2)/2 can written (cf. [1, p. 358]) as 

(2.81) Y v{z) = -!(£) P„(*) + £ log (|) J„(*) - -L (|) " Q v[z) 
where 



e=o 

oo 



(2.83) Q„(z)^!]T(-l) 



£ h e + h v+ i -2-f e (z 



•21 



where hi = J2k=i (^o = 0) and 7 e is the Euler constant. Thus, the kernel 
involves two types of singularities for n > 4 and can be factored as 

(2-84) K k (r) = ^+p k (r)log(r) + K k (r) 

where 

(2.85) a k (r) = ^S, 



(2.86) p k (r)=p k (0)A n (kr) where ^ {Q) = _ _^_ 



n-2 



and 

(2-87) K k (r) = , ' ... ( * V -hog ( k - ) \. lA kr) + , 



4r(n/2) V2v^ 



tA„(fcr)J. 
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The limiting values of and are given by 
(2.88) a k n {Q) = ^fS, a k (0)=0 

(2-89) # +2 (0) = ^H£(0), /? 2 fe (0)--^ 

The following table summarizes formulae and values of and /3 k for the first few 
even dimensions. 



n 




a*(0) 




#(0) 


2 








-J (kr)/(2n) 


-1/(2tt) 


4 


1/(4tt 2 ) 


1/(4tt 2 ) 


-fcJi(fcr)/(47r 2 r) 


-fc 2 /(8^ 2 ) 


6 


(l + (W2) 2 )/(4^ 3 ) 


1/(4tt 3 ) 


-fc 2 J 2 (fcr)/(87rV) 


-fc 4 /(647T 3 ) 


8 


(2 + (fcr/2) 2 + (fcr/2) 4 /2)/(47r 4 ) 


1/(2tt 4 ) 


-fc 3 J 3 (fcr)/(167r 4 r 3 ) 


-k 6 /(768ir 4 ) 



The value of K k at the origin is given by 

Hence, if*(0) for any even n > 2 can be generated by the recurrence relation 
(2.91) ^(O.—AKoj + j^^^) 

from the initial value 

The existence of the additional logarithmic singularity does not add any diffi- 
culty to the construction of the quadrature weights; now, the formula (2.24) simply 
contains one more term. 



(2.93) w e 



' a k (0) + (3 k (0) ^ + Kl (0) _ I = 

v a*(ry) - 0°(r,)) <p(rt) + P k (n) {f t - <^(r t )) <p{rt) t + 

where <fi a (r) = l/r"~ 2 and <j>^(r) — log(r). The regularized singularities <j>" and <j>^ 
are obtained by the same procedure, independently of each other. 

3. Convergence Analysis 

In this section, we present the rate of convergence of the presented quadrature 
rules depending on the regularity of the data. The main result is that our corrected 
trapezoidal rules applied to smooth data converge faster than any algebraic order of 
accuracy. We begin with defining useful regularity classes. The Fourier coefficients 
of functions in each class shares a common form of upper bounds. Then, we present 
the accuracy of the (uncorrected) trapezoidal rule for the regular integral (2.14). 
Finally, the accuracy of the corrected trapezoidal rule for the singular integral (2.18) 
is presented. 
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3.1. Decay characteristics of Fourier coefficients. In spectral contexts, the 
classification of functions by the decay characteristics of their spectral coefficients 
will always be the best one since the error of a numerical scheme is in principle 
determined by the rate of decay of the coefficients. However, in common situations, 
an a priori estimate is not likely to be available. On the contrary the smoothness 
of a function is more accessible, and upper bounds for the Fourier coefficients |-F(fc)| 
can be obtained by repeated applications of integration-by-parts. One of the earli- 
est application of this technique (in a somewhat diffrcrcnt direction) can be found 
in [16]. We present the extension of one-dimensional results (cf. [5, 6]) to higher 
dimensions. Although this approach has been widely exercised, the employed regu- 
larity conditions vary depending on authors resulting, sometimes, a less tight error 
bound for the trapezoidal rule. We begin with the following definition. 

Definition 3.1. For any non- negative integer P, we denote by C p er the class of 
periodic (or periodized) C p functions with U as a period, which satisfy the following 
conditions. Let F g C p er . 

(1) For P > and P is even, A (P ^F g C° er . 

(2) For P > and P is odd, djA^W F g C° er for all j = 1, ... , m. 

(3) For P = 0, (a) U consists of a finite number of disjoint sub-domains on 
each of which F is C 1 up to the sub-domain boundary and C 2 with AF in 
L 1 , and (b) the boundary of each sub-domain is of C 1 and each connected 
component of the boundary has two adjacent sub-domains. 

Then, an F g C p er is of the Holder class C p,x . With an arbitrary assignment 
of the orientation, unit normal vectors are well-defined on the sub-domain bound- 
aries, and so is the trace of <9(AL p / 2 J F) / 'dn. The conditions on Cp er regarding 
the piecewise smoothness arc to enable the application of the Green's identities. 
Those piecewise smoothness enable a more precise estimate by considering two 
types of manageable singularities in the (P+2)nd derivatives of a function - L 1 and 
if -1 . Take the function F = max(0, 1 — 4i 2 ) on [—1, 1] for example. F is in Cp er 
and its Fourier coefficients are given by F(k) = (sinc(7rfc/2) — cos(7rA:/2))/(7rfc/2) 2 
(0(|fc|~ 2 )). Less rigorously, we can expect the result without the precise evalu- 
ation; the second derivative F^ consists of two delta functions at ±1/2, hence, 
the Fourier coefficients of F^ are 0(1), which results in the 0(|/c| _2 )-decay of 
F(k). On the contrary the function ^/max(0, 1 — 4i 2 ) is C° but docs not satisfy 
the piecewise regularity conditions (hence, is not in Cp er ). Its Fourier coefficients 
are Ji(nk/2)/(2k), which is of 0(|fc|~ 3//2 ). Wc consider the first example can be 
observed more frequently than the second. The second example illustrates that we 
can enrich the classification by introducing (weaker) classes with half-integer indices 
(L 1 first derivatives). However, in this paper, we follow the virtue of simplicity. 

Lemma 3.2. Let F g C p . Then, 



Proof. (1) For P — 0, we apply Green's second identity on each sub-domain. Let 
r be the union of the sub-domain boundaries with arbitrarily assigned orientation. 
Utilizing A{e- M -y) = -7r 2 ||fc|| 2 e -" fe ^ and the continuity of d(e~ i7vk - y ) / 'dn, wc 



(3.1) 
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obt 



am 



(3.2) 2 m F(k) 



^l|fc|| 2 {/r 


"<9FT 




dn J 



Ju 



(y)e- M -y dS y - / AF(y)e- M - y dy 



where J-] is the jump discontinuity across T. Therefore, |.F(fc)| = O (|jfc||~ 2 ). (2) 
For P = 1, we can apply Green's first identity. By the continuity of the first 
derivatives, the boundary integral disappears. 

(3.3) F(k) = Yl k > j u d ^y) e ~ ivk ' y d y- 

Since djF e Cp er by definition, the integral for each j is O (||fc||~ 2 ). By the Holder 
inequality, | ^,kj\ < y/m\\k\\. (3) For P > 1, apply Green's second identity [P/2\ 
times with the boundary integrals vanishing. Then, 

/ -1 \ LP/2J r 

(3.4) F(k) = [^ m ,) J u ^F{y)e- k -ydy. 

Apply the result of (f ) if P is even, and apply (2) if odd. □ 

The upper bound given by Lemma 3.2 is somewhat conservative; for m > 1, the 
integrals in (3.2) contribute to an additional decay (of possibly fractional order). 
However, in one-dimension, the integral on T becomes point-wise evaluations of the 
jump discontinuity, hence, we cannot expect a faster decay. Also, one should note 
that the smoothness is not the only factor which governs the decay characteristics. 
Consider the function max(0, 1 — 4r 2 ) 2 which is in Cp er . Its Fourier coefficients are 
((12-4(7rfc/2) 2 )sinc(7rfc/2)-12cos(7rfc/2))/(7rfc/2) 4 , hence, are of (9(|fc|- 4 ), which 
is faster than the estimated 0(|fc| -3 ). Although they are in different regularity 
classes, the function max(0, 1 — 4r 2 ) 2 € Cp er exhibits the same rate of convergence 
by the trapezoidal rule as a more smooth function max(0, 1 — 4r 2 ) 3 € C 2 er . 

Lemma 3.3. Let F E C£ er . Then, Y<kei™ \F( k )\ <oo if P > m-2. 

Proof. Rewrite the infinite series by |F(0)| + Y^=i^l\\k\\ 2 =n Let r(n) be 

the number of lattice points on the sphere {||A:|| 2 = n}, which is also known as 
the sum-of-squares function. Then, R(n) = J2n=o r ( n ) 1S tnc uumbcr of lattice 
points in the ball {||fc|| 2 < «}, and the Gauss' circle problem in Z m states that 
R{n) = 0{n m / 2 ) (cf. [18]). Since \F(k)\ = 0(n-( p + 2 )/ 2 ) for ||fc|| 2 = n by Lemma 
3.2, the series is bounded above by F(Q) + CY^ = \ r(n)/n^ p+2 ^ 2 for some C. Since 
r(n)/n( p + 2 )/ 2 = (R(n) - R(n - I))/n( p + 2 )/ 2 = 0{{n m ' 2 - (n - l) m / 2 )/n( p + 2 )/ 2 ) 

and _ ( n _ l)m/2 = Q^m/2-1^ r{n y n (P+2)/2 = ( n -(P-m+4)/2^ which 

implies that the series converges if P > m — 2. □ 
Corollary 3.4. Let F e C p r . The multiple Fourier series 

(3.5) J2 F(k)e t7rk - y 

converges uniformly to F(y) if P > m — 2. 

The criterion P > m — 2 is somewhat excessive in higher dimensions. As men- 
tioned earlier, if we consider the additional decay, the convergence criterion can be 
relaxed. 
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3.2. Accuracy of the trapezoidal rule. By substituting the Fourier series (3.5) 
into the definition of DFT (2.4), we obtain 



(3.6) F k = -±=Y1 F{k')Y.^ (k '' kyVl - 

The inner sum is 2 m N if k'j = kj + 2qjNj for all j and for some q £ Z m , and 
vanishes otherwise. Therefore, 

= ^( fcl + 2qiNl > ■ ■ ■ ' km + 2 1raN m ) 

(3.7) =F(k)+ J2 F(ki+2q 1 N 1 ,...,k m + 2q m N m ). 

9 ez™\{o} 

The above equation describes that the coefficients Fk obtained from the DFT con- 
tain the aliasing error from the harmonics (X^o)- 

Notice that Fq is the integral of F on U obtained by the trapezoidal rule (for 
periodic functions, hence, considering only one side of the boundary), and F(0) is 
the exact integral. Hence, an error bound of the trapezoidal rule is given by 

(3.8) \F(0)-F \< \F(2q 1 N 1 ,...,2q m N m )\. 

9 ez™\{o} 



Let N m i n = min(iVi, . . . , N m ). The following theorem states the order of accuracy 

<p 

per 



of the trapezoidal rule applied to a C p function. 



Theorem 3.5 (The trapezoidal rule). For any F £ C£, r with P > m — 2. The 



J per 

error of the trapezoidal rule applied to F is 0(N~^ +2 ' > ). 
Proof. Let p = {2q 1 N 1 , 2q m N m ). Then, ||p|| > 2JV min ||g|| and 

\F(P)\ < C/(N min \\q\\y p +V 
for some C by Lemma 3.2. Rewriting J2 qe z™\{o} h V EiT=i E|| 9 p=n' 

A, rain n —l 

where r{n) is defined in the proof of Lemma 3.3 and is 0(n" 1 / 2-1 ). The series on 
the right side converges if P > m — 2. □ 

3.3. Accuracy of the singular trapezoidal rule. From (3.5) and (3.7), the 
error of the interpolation F defined by (2.7) can be written as 

F(y)-F(y)= ]T F(k) e™^ 
(3-9) ^ _ . 
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Hence, the error of the singular quadrature rule (2.18) is given by 

< e \F(k)\\m\ 



(3.10) 



/ (F(y)-F(y))0(r)dy 
Ju 



k0 N 

+ E ^ (fc)| S + 2<ziiVi ' • • ■ > km + 2 ^ N ^\ 

keiN gez m \{0} 

Theorem 3.6 (The singular quadrature rule). Let F e C p er with P > m — 2. TTie 
error o/ the quadrature rule (2.18) is 0(N^ n < -^ +2 ~^ m ~ 1 ^ 2 ^). 

Proof. (1) First, consider the first summation in (3.10). Since F(k) — (3(||fc||~( p+2 )) 
by Lemma 3.2 and <j>(k) = 0(\\k\\-^ n+1 ^ 2 ) by Lemma 2.1, 

E \ F ( k )W&( k )\ -EE n (P+m/2+5/2)/2 = X! n (P+m/2+5/2)/2 
fc^liv n=iVV„ ||fc|| 2 =n ™ =JV min 

where r(n) = 0(n m / 2_1 ) as in the proof of Lemma 3.3. Thus, the right side is 
bounded by C£n=N 2 n -^ p - m / 2+9 ^/ 2 , which converges if P > (m - 5)/2 and 

is °(^mm +:Mm ~ 1)/2) )- ( 2 ) Consider the second term in (3.10). Since \kj\ < Nj 
for j = 1, . . . , ra, the minimum of the parabola (kj + 2qjNj) 2 is min((2qj ± 1) 2 N 2 ) 
if qj ^ 0. Hence, for qj ^ 0, (kj + 2q J N J ) 2 > (2\qj\ - 1) 2 N 2 > q 2 N 2 , and the last 
inequality holds for qj = also. Therefore, ||(fci + 2q\N\, . . . , k m + 2q m N m )\\ > 
\\q\\N min and \F(h + 2q 1 N 1 , . . . , k m + 2q m N m )\ = 0((N min \\q\\)-^ p + 2 )) by Lemma 
3.2. Then, the second term is bounded by 

(3.11) c\\m\+f: E «- (m+1)/4 ){^L p+2) £ E p- {p+2)/2 }> 

^ n=l \\kP=n ' p=l || 9 || 2 =p ' 

where we utilized Lemma 2.1 to get the upper bound for |(/>(fc)|. The infinite series 
in the second pair of braces converges if P > m — 2. Since 

(3.12) E E ^ (m+1)/4 < C E n (m - 5)/4 - 0(N^- 1)/2 ), 

n=\ ||fe||2=n n=l 

the second term in (3.10) is 0(N m { ^ +2 ~ {m ~ 1)/2) ). Thus, (1) and (2) complete the 
proof. When the kernel is less singular and \<ft(k)\ decays faster than (3(||fc||~( m+1 )/ 2 ), 
(3.11) becomes dominant and determines the rate of convergence. If \4>(k)\ decays 
faster than 0(||fc|| _m ), the first term in (3.11) is 0(1) and the rate of convergence 
is the same as the regular trapezoidal rule, i.e. (P + 2). □ 

Combined with Theorem 3.5, Theorem 3.6 indicates that the order of accuracy of 
the corrected trapezoidal rule presented in this paper is at least (P + 2 — (in — 1) /2) 
for the integral (1.1) with / e C p er . If / is smooth, the corrected trapezoidal rules 
converges faster than any algebraic order. The dimension-dependent degradation 
from the intended (P + 2) (as the regular trapezoidal rule) by the amount of — (m — 
l)/2 is more like a technical outcome; in actual experiments with m < 3, we have not 
experienced any obvious degradation in the rate of convergence and results strongly 
imply that the corrected trapezoidal rule exhibits the same rate of convergence 
as the usual trapezoidal rule without the singular kernel. The same is true for 
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the condition P > m — 2. Our cautious conjecture is that, as briefly mentioned 
previously, there is an additional decay for m > 1 originating from the integrals in 
(3.2), which cancels (at least, a part of) the degradation. 

The intention of the somewhat complicated proof of Theorem 3.6 can be well 
illustrated by considering simple one-dimensional cases. For C^ er functions, both 
the trapezoidal rule and the corrected trapezoidal rule exhibit [P + 2)nd order con- 
vergence. The interpolation error, however, is only of (P + l)st order. The mutual 
cancellation of overshooting and undershooting, which has been well explained for 
the regular trapezoidal rule, holds identically for our singular quadrature rules. In 
the singular cases, the additional one in the rate of convergence originates from the 
decay characteristics of |0(fc)|. 

4. Numerical Examples 

4.1. Test of convergence: Helmholtz kernels with k — 0. We evaluate the 
convolution with the non-oscillatory kernels on R m (m = 1, 2, 3): 

(4.1) u = K°*f 

for n = m and m+ 1. For n = m, u corresponds to the volume potential induced by 
the surce /. For n = m + 1, the convolution corresponds to the application of the 
single layer operator as a boundary to boundary integral operator on m-dimensional 
flat boundary in M m+1 . Convolutions were performed on uniform grids with the 
spacing 6/iV on the physical domain [—3, 3] m . The errors are reported in norm 
from the values of the numerical and the exact solutions obtained at grid points. 

Table 1 and 2 summarize the results of one-dimensional convolution with the 
logarithmic kernel: = — log(r)/(27r). Reference solutions are computed up 
to the machine precision using adaptive quadratures. For the demonstration, we 
selected three different sources, 

/ G (r)=cxp(-(r/a) 2 ) (a = 1/2) 

(4.2) f B (r) =exp (12 - 12/(1 -(r/a) 2 )) (a = 2) 

f P (r)= max (0, 1 - (r/a) 2 ) 7 (a = 2), 

based on their regularity characteristics. The Gaussian (fg) is analytic but is 
not compactly supported. However, on the boundary of the domain employed, 
I/gI ~ 10~ 16 and any significant error of the domain truncation has not been 
observed. Actually, the spectrum of the Gaussian (in any dimension) exhibits the 
most rapid decay, which makes it a perfect specimen for the test of the spectral 
accuracy. Table 1 shows a clear super-algebraic convergence of our quadrature 
rule; each refinement doubles the number of correct digits. The bump function 
Jb is smooth and compactly supported, hence, satisfies our assumptions faithfully. 
However, its Fourier coefficients decays quite slowly compared with the Gaussian. 
The left column of Table 2 shows that the estimated order keeps increasing but the 
increase is slower than that of the Gaussian cases in Table 1. The function fp is 
not smooth but in Cp er . Hence, the estimated order is fixed to the algebraic order 
of 8 (right column in Table 2), which is expected in the convergence analysis. 

Recall §2.4 that we can use arbitrarily refined grid during the construction of 
the quadrature weights. Such a refined construction increases only the construction 
time but does not affect the convolution time since the obtained weights will be 
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fa (without ref.) fa (with rcf.) 



N 


En 






En 




l°S2^f 


5 


5.58 x 10" 


2 




5.59 x 10" 


-2 




10 


3.26 x 10" 


3 


4.1 


3.26 x 10" 


3 


4.1 


20 


1.30 x 10- 


6 


11.3 


1.30 x 10" 


6 


11.3 


40 


3.32 x 10" 


13 


21.9 


3.89 x 10" 


■16 


31.6 



Table 1. Loa errors of (K° * fa) on R: (left) without and (right) 
with using the doubly-refined grid for the construction of the 
weights. 



truncated in the frequency domain. The objective of the use of a refined construc- 
tion grid is to suppress errors originating from insufficient samplings of a, ip, and 
K. Since K® are non-oscillatory (a is a constant and K = 0), the error during the 
construction originates mostly from ip. In Table 1, the results on the left column 
are obtained using the weights constructed on the same grid as the data. For the 
results on the right column, we doubled the sampling frequency of the construc- 
tion grid. For small N, the difference is negligible since the error from the data is 
dominant. However, the interpolation error of the Gaussian decreases more rapidly 
than that of p. Hence, if we use the same grid, the construction error stand out 
eventually. The remedy is simple; we can use a slightly higher sampling frequency 
for the construction. As for fg and fp, their interpolation errors decrease more 
slowly than ip, and the results are almost indistinguishable with or without using 
a refined construction grid (see Table 2). 

Table 3 shows the results for (K® * fa) m higher dimensions (to = 2,3 and 
n = m, m + 1). For the Gaussian, we know the exact solutions, which are given by 
(4.3) 



(to = 2,n = 2) u(r) = 

(to = 2, n = 3) u(r) = 

(to = 3, n = 3) u(r) = 

(to = 3, n = 4) u(r) = 




where p — r/a and a = 1/2. As in the one-dimensional case, we can observe clear 
spectral rates of convergences. 

4.2. Test of convergence: Helmholtz kernels with k ^ 0. Convergence tests 
were performed for the Helmholtz kernels K% with nonzero wavenumber. First, we 
present results with real wavenumber k — 2n on [— 3, 3] m . The Gaussian Jg was 
selected as the source. Since the Fourier coefficients of fa decay rapidly, the quad- 
rature rule, if is is properly constructed, should exhibit a similar fast convergence. 
Unlike non-oscillatory kernels of the previous section, we do not know the exact 
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!b fp 

N E N log 2 %f E N log 2 



5 


4.17 x 10" 


2 




1.46 x 10" 


-2 




10 


7.21 x 10" 


-4 


5.8 


5.65 x 10" 


■5 


8.0 


20 


1.45 x 10" 


■6 


9.0 


2.36 x 10- 


-7 


7.9 


40 


9.25 x 10" 


10 


10.6 


7.31 x 10" 


-10 


8.3 


80 


2.36 x 10- 


-14 


15.3 


4.33 x 10" 


-12 


7.4 



Table 2. L x errors of (K$ * f B ) and (K$ * f P ) on R. 



N 



(to = 2,n = 2) 



E 



N 



5 1.06 x 10- 1 

10 3.96 x 10~ 3 

20 8.99 x 10~ 7 

40 5.55 x 10~ 16 



log. 



E 



N/2 

2 



4.7 
12.1 
30.6 



(m = 2,n = 3) 



AT 



4.88 x 10~ 2 
4.70 x 10" 3 
2.35 x lO" 6 
3.33 x 10" 16 



log 



2 gjv 

3.4 
11.0 
32.7 



N 


(to = 


3,n 


= 3) 


(m = 


3,n 


= 4) 












log 2 £ # /2 


5 


4.04 x 10- 


2 




1.41 x 10- 


2 




10 


4.10 x 10- 


3 


3.3 


5.03 x 10- 


3 


1.5 


20 


1.19 x 10- 


6 


11.8 


3.22 x 10- 


6 


10.6 


40 


1.05 x 10- 


15 


30.1 


3.05 x 10- 


16 


33.3 



Table 3. errors of (K° * fa) on R' 



solution even for the Gaussian. Moreover, it is very difficult and expensive with 
adaptive quadratures to obtain the reference solutions accurate up to the machine- 
precision. Hence, the reference solutions were evaluated only at the origin, where 
the error is likely to be largest. The results are summarized in Table 4 and 5, from 
both of which we can observe the expected spectral rates of convergences. 

Similarly to Table 1, Table 4 illustrates the enhancement in the accuracy by 
using a higher sampling frequency for the construction grid. However, the results 
with the oscillating kernels exhibit quite different aspects: first, the level of the 
construction error is significantly larger than the previous non-oscillatory cases. 
Second, the corruption occurs also for small N. Recall, for non-oscillatory kernels 
in the previous section, the main source of the construction error was ip, whose 
Fourier coefficients decay quite rapidly. Although the interpolation error of the 
Gaussian exhibits faster decay than that of tp, the Gaussian is a somewhat special 
case. For many other functions (such as fg and fp), ip is good enough not to cause 
such a issue. However, for oscillatory kernels, the major source of the construction 



24 



JAE-SEOK HUH AND GEORGE FANN 



(to = 1, n = 2, without ref.) (m = 1, n = 2, with ref.) 



5 


4.11 x 10" 


2 




7.95 x 10" 


2 




10 


4.66 x 10" 


2 


-0.2 


6.47 x 10" 


3 


3.6 


20 


2.89 x 10" 


-4 


7.3 


2.82 x 10- 


■6 


11.2 


40 


2.61 x 10" 


-11 


23.4 


3.93 x 10" 


■17 


36.1 



Table 4. L M errors of * fa) on R with k = 2n: (left) without 
and (right) with using the doubly-refined grid for the construction 
of the weights. 



(to = 2,n = 2) 



(to = 2,n = 3) 



iV 



AT 



5 4.01 x 10~ 2 

10 1.14 xlO" 2 

20 2.46 x 10~ 6 

40 2.08 x 10 



17 



1.8 
12.2 
36.8 



1.02 x lO" 1 
1.26 x lO" 2 
4.77 x 10~ 6 
2.55 x 10~ 16 



3.0 
11.4 
34.1 



(to = 3, n = 3) (to = 3,ti = 4) 

N E N log 2 %^ E N log 2 



5 


4.61 x 10" 


2 




1.13 x lO" 1 




10 


1.52 x 10" 


2 


1.6 


1.81 x lO" 2 


2.6 


20 


2.95 x 10" 


6 


12.3 


6.17 x 10~ 6 


11.5 


40 


2.96 x 10- 


■17 


36.5 


4.13 x 10~ 16 


33.8 



Table 5. errors of (K* * f G ) on R m with k = 2rr. 



error is a and K. Simply, we need a construction grid fine enough to represent the 
oscillating a with similar accuracy as the data (K oscillates similarly to a). Hence, 
the sampling frequency of the construction grid should be increased proportionally 
to the wavenumber. 

On the other hand, if the sampling frequency of the data is not high enough 
for the representation of the oscillating kernel, the sampling frequency is not high 
enough for the representation of the solution either. Hence, ironically, in most 
of practical applications such as examples in §4.3 and §4.4, we cannot expect the 
enhanced accuracy by merely increasing the sampling frequency of the construction 
grid. The first convolution from the exactly given data may be accurate at each grid 
point. However, the interpolation error of the first solution will not be as accurate 
as that of the data, due to the insufficient grid resolution. This first solution with 
insufficient accuracy will be (a part of ) the source of the next convolution. Hence, 
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Im k 

Figure 1. The Helmholtz kernel K\ on E: (left) the weights for 
k — 4d in the frequency domain and (right) errors up to k = 50i. 

if an appropriate grid is chosen for the data and the solution, we can use the same 
grid for the construction. 

The situation becomes more complicated if Im(fc) 7^ 0. Take K\ with a pure 
imaginary wavenumber k = iX for example. 

(4-4) Kl(r) = X>(iAr) = *?M „ Xi"-* 

4 27T V o7rAr 

As described in §2.6, the kernel can be factored as K\{r) = /3| (r) log(r) + K%(r) 
with 

«,5, flJW ._*^)._^„^ ( ^_ fa - ) . 

Jo and -Ko are the modified Bessel functions of the first and the second kinds. 
Notice that the exponentially decaying kernel is decomposed by two exponentially 
increasing fi\ and K\. With finite precision arithmetic, the meaningful signal 
(~ e~ Xr /*Jr) in 0% IS completely lost for Xr ^> 1. Thus, we should choose sufficiently 
small Br so that the /3f evaluated within Br remains small ~ 1. The left figure in 
Figure 1 shows the spectrums of two sets of the weights generated for A = 4. For 
the accurate weights, the radius of Br is set to 1. The inaccurate weights obtained 
with Br of the radius 3 deviate wildly after the mid-frequency. As we increase A 
with the fixed radius of 3, the points of deviation move to lower frequencies, and 
soon only the noise remains. Since the construction grid should be fine enough 
to resolve ip with Br, its sampling frequency should be increased accordingly. In 
order to prevent the generation of a huge construction grid due to the high sampling 
frequency, we can take only a small subset of grid containing Br. Outside of the 
subset, the kernel is practically zero, hence, the spectrum of the weights can be 
obtained by zero-padded DFT followed by the truncation of the unnecessary high 
frequency results. 

In the right figure of Figure 1, the errors with the appropriately constructed 
weights are less than 10~ 14 . Without the A-dependent control on Br, the error 
rapidly increases after A > 5. 
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N 


En 




l°S2^f 


80 


1.42 x 10- 


-l 




160 


2.08 x 10- 


4 


9.4 


320 


2.07 x 10- 


-7 


10.0 


640 


7.42 x 10- 


■11 


11.4 



Table 6. -Loo error of the solution of (4.8): k — 5ir on [— 6,6] 2 . 



4.3. Application: Lippmann-Schwinger equation. Consider the acoustic scat- 
tering problem in an inhomogeneous medium. Let the constant k be the wavenum- 
ber of the medium at infinity (or we may call it the ambient wavenumber) and let n 
the non-constant refractive index such that (n — 1) is compactly supported. Then, 
for time-harmonic problems, the acoustic pressure u satisfies the equation, 

(4.6) Au + k 2 n{x) u = 0. 

Let u l be the given incident wave which satisfies Au 1 + k 2 u l = 0. We assume 
the scattered field u s = u — u L satisfies the radiation condition at infinity. The 
explicit form of the radiation condition depends on the dimension. In this paper, 
we consider an example in R 2 , where the radiating scattered field satisfies 

- / du s \ 

(4.7) lim y/r iku = 0. 



dr 

The above scattering problem is equivalent to solving the integral equation, 

(4.8) u{x) - k 2 f K*(r) (n(y) - 1) u(y) dy = u\x), 

JD 

which is also known as the Lippmann-Schwinger equation. The linear operator on 
the left side of the equation (denoted by C) can be written as 

(4.9) C = I-k 2 KM 

in terms of the convolution JCu = K\ * u and the multiplication Afu — (n — l)u. 
Since /C is compact and M is bounded, £ is a Fredholm operator of the second kind 
and the problem is well-posed (cf. [9, §8]). 

In this example, we consider a medium with three bumps, 

(4.10) n(y) = 1 - 0.9 ]T ^H^vi-aO'-d*-*) 3 )- 1 ) , 

with the centers at (a,,6j) = (1,0), (—1,3), and (—1,-3). Inside of the bumps, 
the wavespeed slows down. The equation Cu = u l is solved using GMRES for 
incident planewave with the wavenumber k = 5n and the direction (1,0). The 
domain of the computation is [—6, 6] 2 . In order to measure the error, the reference 
solution is computed on the uniform 1280 x 1280 grid. For the convergence test 
presented in Table 6, we measured error on the grid points. Note that, unlike 
the convergence tests in previous sections, the reported errors are not just for the 
convolution itself but includes the combined effect of the whole solution procedure. 
Figure 2 depicts the total field the scattering problem. 
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4.4. Application: boundary integral equation. The presented quadrature rules 
can be applied for integral operators on curves in M m . Let each curve be repre- 
sented by a smooth periodic parameterization: s £ [0,L] H> y(s) € K m . Assume 
y'(s) 7^ for every s G [0,L]. For s,i e [0,L], define the ratio, 

(4.11) 7(M) = ^| where r(s,t) = \\y(t) - y(s)\\. 

\t-s\ 

Then, 7(i,t) = ||y'(i)|| and j(s, •) is smooth, positive, and periodic on [Q,L]. The 
logarithmic and the power-law singularities can be rewritten as 

(4.12) log(r)=log(|«-t|)+log(7) 

(4.13) r~ v = ~f- v \s - t\~ v . 

Thus, the above radial singularities in K m can be recast to singularities of the same 
kind on [0,L]. 

In this example, we solve the exterior scattering problem, 

(4.14) l - i>(s) + J D(s, t) m \\y'(t)\\ dt - ikJ^S{s, t) ^(t) \\y'(t)\\ dt = -«*(«), 

where the combine integral equation is employed to avoid the interior resonance 
issue. The single layer kernel S is identical to the Hclmholtz kernel K\ except the 



28 



JAE-SEOK HUH AND GEORGE FANN 



additional log(7) in the smooth remainder; 

l -H^ (kr) = a s (s, t) logfla - t\) + S(s, t) 

Jo(kr) 
2ir 

1 

i leuler 1 , / k\\y' (t)\\ 

4 2tt 2tt ° g ^ 2 
where jeuier is the Eulcr constant. The double layer kernel D is given by 

_ d fi A [ "(*) • (V(t) ~ y{s)) 



(4.15) 


S(s,t) 


(4.16) 


a s (s,t) 


(4.17) 


a s (t,t) 


(4.18) 


S(t,t) 



D(s,t) = n(t) ■ Vtf£(r) = - ^-H™(kr) } 
(4.19) = k ^{Y 1 (kr)- l J 1 {kr^^ {t) -^ t) - y{s)) 



4 

where n(t) is the outward normal vector at t. Since (cf. §2.6.2) 
(4.20) zYi(z) = 2 -^l\o g {z) - - + 0(z 2 ), 

7T 7T 

the kernel can be factored by 



(4.21) 


D{s,t) 


= a d (s,t) lo{ 


5(1*- 


t|) + £)(«, t) 




(4.22) 


a d (s,t) 


krJi{kr) 
2tt 


fn(t) 


■(»(*) 

r 2 j 


> 


(4.23) 


a d (t,t) 


= 








(4.24) 


D(t,t) 


= lim i 

2ir s-n \ 


'»(*) 


■(y(*)-y(«))l 


c(t) 




r 2 J 


2^ll2/'WII 2 



where c(t) is the curvature at t. The above information is all we need to construct 
the quadrature weights (at each target point). Beware that a fast convolution 
cannot be applied for this case since the kernels are not functions of (t — s), hence, 
the resulting discrete operator is not a circular matrix. 

Figure 3 shows the total field constructed from the obtained solution i[>. Each 
scatterer is a translation of the popular kite shape, 

(4.25) y(t) = (cos t + 0.65 cos(2f) -0.65, 1.5 sin t) ie[0,2?r]. 

The incident planewave with k — 5ir comes from the upper-left corner of the domain 
[—8, 8] 2 to the lower-right corner. For the interaction between disconnected curves, 
the kernels become smooth periodic, hence, the usual trapezoidal rule can be used. 

Table 7 shows the error of the solution ip on the boundary and the error of the 
far-field evaluate from the obtained tp. The far-field signature can be computed by 



—iir/4 r 

(4.26) Uoo (x) = ^= / {kn{t)-x + k)e-^-y {t ^{t)\\y'{t)\\dt 

VSirk Jr 



-iir/A 

It 

where x — x 1 1 x 1 5 X G M. 2 . The reference solution is obtained with N — 640. 
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solution error 



far-field error 



N 



E 



N 



En 



80 1.16 x 1CT 1 
160 1.61 x 10" 5 12.8 
320 4.17 x 10~ 14 28.5 



1.45 x lO" 1 
2.06 x 10~ 7 
1.34 x 10~ 14 



19.4 
23.9 



Table 7. errors of the exterior scattering problem (4.14): k = 5ir. 




Figure 3. Total field of the exterior scattering problem (4.14) on 
[—8, 8] 2 : k = 5tt. The obstacles are centered at ±2 on x axis. 



4.5. Conclusions. We have presented the construction scheme for corrected trape- 
zoidal rules for integral operators with weakly singular kernels. Numerical results 
show that the quadrature rules converge as fast as the trapezoidal rule applied for 
non-singular periodic integrands as predicted by the presented convergence analysis. 
Especially, for smooth data, the quadrature rules exhibit super-algebraic conver- 
gence. The examples in §4.3 and §4.4 demonstrate the effectiveness of the proposed 
quadratures in the high precision solve of integral equations. 

Appendix A. Asymptotic Expansion of 
Dividing the domain of integral, we rewrite M± (p) as 

1 

(A.l) (p) = (a) + Up, p) where Up,, p) = — I t^ 1 cos(f) dt. 

P M J a 
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Applying intcgration-by-parts (27V) times, we obtain 

(A.2) Ia{fljp) = I N { ^ p) + t3l^ {e _ fl) [ P t»-i-™ cos{t )dt, 

P «_i Jo 



where 

I?(p,p) =p- 1 Mp)P N (P,P)-a- 1 (*/pYsm(a)P N (^a) 
+ p- 2 cos(p) Q N (p, p) - a- 2 (a/pY cos(a) Q N {p, a). 
Polynomials P N and Q N are denned by 

JV-l 21 

(A.4) P N (p,t) = Cit~ u where C e = (-l) e " m) 

(A.5) Q N (p,t) = (p-l)P N (p-l,t). 

Then, the relative error is given by 

(A.6) \l a ( Ji ,p)-I?( Ji ,p)\<^\I a (p)\. 

By choosing a = lAn and N = 13, the relative error is less than 10~ 16 . Moreover, 
with a being an integer multiple of n, l£ has a slightly simpler form. 

Appendix B. Asymptotic Expansion of 
Similarly to , we divide the domain of integral; 

(B.l) M^(p) = M^(a)+I a (p,p) where I a (p, p) = \f V 1 J (t) dt 

P J a 

Utilizing t Jo (t) = (tJ\{t))' and J\{t) = —J (t), (2N) applications of the integration- 
by-parts result in 

(B.2) I a (p,p) = I?(p,p)+ { —^- \\{2(.-p) 2 / t^- 2N J (t)dt 

where 
(B.3) 
with 

iV-l l 

(b.4) p N {p,t) = c * r2£ whcre °i = - ^ 

e=o j=i 

N-l 

(B.5) Q N (n,t)= Ysip-^-^Cit- 21 . 

1=0 

The relative error of the asymptotic expansion is given by 
(B.6) \l a (p, p) p)\ < -^2- \Up)\ ■ 

With a = 44.7593189976528217 (which is the 14th zero of J x ) and N = 15, the 
relative error is less than 10~ 16 . 



I?(p,p)=p- 1 J 1 (p)P N (p,p)-a- 1 (a/prj 1 (a)P N (p,a) 
+ p- 2 Mp) Q N (P, P) ~ a- 2 (a/p)» J (a) Q N (p, a) 
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